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Abstract 

Two Galerkin methods are applied to a problem of convection with uni- 
form internal heat source are given. With each method analytical results are 
obtained and discussed. They concern the parameter representing the heating 
rate. Numerical results are also given and they agree well with the existing 
ones. 
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1 The physical problem 

Natural convection induced by an internal heat source is a phenomenon which 
has been intensively studied, especially in order to point out its influence on other 
processes. The motion in the atmosphere or mantle convection are two among such 
phenomena [15] . They bifurcate from the conduction state as a result of its loss of 
stability. A major importance is given to thermal convection processes in terrestrial 
bodies driven by internal heat sources in which the heat source is a function of 
time and, moreover, can vary from one terrestrial body to another. In spite of 
their importance, due to the occurrence of variable coefficients in the nonlinear 
partial differential equations governing the evolution of the perturbations around 
the basic equilibrium, so far these phenomena were treated mostly numerically and 
experimentally. In \A\ we carried out a linear study for the eigenvalue problem 
associated with the equations for a convection problem with an uniform internal heat 
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source in a horizontal fluid layer bounded by two rigid walls [15]. Our method was 
based on Fourier series expansions for the unknown functions. Numerical results and 
graphs were given showing a destabilizing effect of the presence of the heat source. 
In [5J another two methods based on Fourier series expansions (a Chandrasekhar 
functions - based method and a shifted Legendre polynomials - based method) were 
used to study analytically the eigenvalue problem deduced in [I]. 

In [T3] a linear stability analysis for a natural convection problem induced by 
internal heating is performed in order to point out the effects of the heat distribution. 
This is a function of both the critical Rayleigh number and the critical wavenumber. 
Some non-uniform distributions were considered along with the uniform one. It was 
shown that a concentration of the heat source near the bottom boundary implies 
a decreasing of the stability domain; namely it lowers the temperature difference 
at which the convection sets in. The variation of the critical wavenumber is small 
and there is only a slight influence of this distribution on the size of the convection 
cells. When the heat source is placed near the top boundary an enlargement of the 
domain of stability occurs. 

Another analytical study for a problem of convection in a fluid saturated porous 
layer heated internally and in the presence of a linearly varying gravity field is 
presented in [8]. It was proved that the principle of exchange of stabilities holds 
as long as the gravity field and the integral of the heat source have the same sign. 
Convection in a medium with internal heat source was also analyzed in [I] by linear 
stability methods and nonlinear stability (energy type) methods. Numerical bounds 
for the critical value of the control parameter, the Rayleigh number, were given and 
the continuous dependence of the solution of the initial boundary value problem on 
the internal heat source was proved. 

In [15] a horizontal layer of viscous incompressible fluid with constant viscosity 
and thermal conductivity coefficients is considered. The performed numerical inves- 
tigation concerned the vertical distribution of the total fluxes and their individual 
components for small and moderate supercritical Rayleigh number in the presence of 
a uniform heat source. In this context, the heat and hydrostatic transfer equations 
are [15] 

d PB 

= ~ PB9 ' (2) 
where rj = const, is the heating rate, 8b, Pb and ps are the potential temperature, 
pressure and density in the basic state. In the fluid, the temperature at all point 
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varies at the same rate as the boundary temperature, so the problem is characterized 
by a constant potential temperature difference between the lower and the upper 
boundaries A9b = 9b — Qb x - Taking into account ([1]) this leads to the following 
formula for the potential temperature distribution [15] 



A9 B / h\ rj 
eB = eB °-—\ Z+ 2) + 2k 



Z ~>2 



(3) 



In nondimensional variables the governing system of equations is 



( dll 

= -Vp' + AU + Gr9% 

dt 

divU = 0, (4) 
dff 

— = (1 - Nz)Uk + Pr~ l A9\ 

where U = (u, y, w) is the velocity, 9' and p' are the temperature and pressure 
deviations from the basic state [6], Gr is the Grashof number, Pr is the Prandtl 
number and N is a dimensionless parameter characterizing the heating (cooling) 
rate of the layer. 

The boundaries are considered rigid and ideal heat conducting, so the boundary 
conditions read 

1 l 

U = 0' = o at z = — and z = -. (5) 

In jl] in order to deduce the eigenvalue problem we considered the viscous in- 
compressible fluid confined into a rectangular box bounded by two rigid walls: 

1 1 

V : < x < di, < y < a,2, — < z < -. We assumed that any unknown 
function in (HI) is of the form from [7] 



f(x, y, z) = F(z)exp[ i[ 2-Km' h 27m — 

V V ai do 



a\ a 2 

2nm' 2txv! , L I r > 

m = ,n = , where a\ = — , a 2 = — , L and / are the box sizes. Here 

ax a 2 H H 

vnl > 1 and n' > 1 are the number of cells in the x and the y direction. 

Another possibility is to assume disturbances periodic in x (period 2n/a) and y 
(period 2%/P), with a growth rate a, also of the form 



f(x, y } z) = F(z)exp[ at + iax + i/3y 
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In this subsequent investigation will concern the condition in which the 

principle of exchange of stabilities is valid. 

In this paper, we treat only the stationary case and this implies that the principle 
of exchange of stabilities is valid. We complete our analytical study from [I], [5] 
with some remarks on the spectral methods used to solve the eigenvalue problem 
governing the linear stability of the basic state for the convection problem with 
uniform internal heat source. 

The eigenvalue problem associated with the equations flU)-© in a horizontal 
fluid layer bounded by two rigid walls, governing the stability of the basic motion 
against normal mode perturbations, deduced by us in [I] has the form 

rp 2 -a 2 )W = e, 

\ (D 2 - a 2 )0 = -a 2 R{l - Nx)W K ) 

with the boundary conditions 

W = DW = = at x = ±~. (7) 

Here the Rayleigh number R = Gr ■ Pr represents the eigenvalue, while W, 0, the 
amplitudes of the perturbations for the velocity and the temperature field respec- 
tively, form the corresponding eigenvector is (W, 0). 



2 On the convergence of the Galerkin method 

In this section, we reveal some aspects of the convergence of the Galerkin method, 
one of the most used method for converting a differential operator boundary value 
problem to a discrete one. 

There are more than one analytical possibilities to solve the system 0H])-([7j). 
However, some remarks on the convergence of the system are in order. First, let us 

perform a translation of variables z = x H — , such that the problem ([6]) becomes 

/ (D 2 — a 2 ) 2 W — = 0, 

X {D 2 — a 2 )0 + a 2 R(Ni — Nz)W = 0, 1 J 

N 

with Ni = 1 + — and the boundary conditions 

W = DW = = at z = and 1. (9) 
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The equations from (jSJ) can be considered as a particular case of a more general 
eigenvalue problem with variable coefficients j3] 

/ (D 2 -a 2 ) 2 W = f(z)&, 

\ (D 2 -a 2 )Q = -a 2 Rg(z)W, [W) 

on < z < 1. 

The mathematical problem reads: for given f(z) and g(z) (in our case f(z) = 1 
and g(z) = N\ — Nz) determine the minimum real positive R over all real positive 
a for which there exists a nonnul solution of the system / fiOj) -/fPj). 

Following Kolomy [9] the convergence of the Galerkin method can be consid- 
ered for the sixth-order equation (D 2 — a 2 ) 3 W = —a 2 R(Ni — Nz)W obtained by 
eliminating G between the two equations from (ITUj) . The following result holds: 

Proposition 1. The operator L = (D 2 — a 2 ) 3 is not symmetric in the sense of an 
L 2 (0, 1) inner product on a space of functions satisfyingW = D 2 W = (D 2 —a 2 ) 2 W = 
at z = 0,1. 

In order to prove Proposition 1, consider the inner product (LW, W*) in L 2 (0, 1) 
with W, W* functions from VL, 

VL := {U G L 2 (0, l)\U = D 2 U = (D 2 - a 2 ) 2 U = at z = 0, 1}. 

The operator L is said to be symmetric if (LW, W*) = (W, LW*) for any W, W* E 
VL. In our case, by direct integration by parts it can be proven that (LW, W*) = 
(W,LW*). However, W* is not a function from VL, namely W* satisfies boundary 
conditions of the type 

W* = D 2 W* = D(D 2 - a 2 )W* = at z = 0, 1, (11) 

whence Proposition 1. In [1] the quoted sixth order equation together with the 
boundary conditions was investigated using spectral methods based on trigono- 
metric Fourier series and good numerical results were obtained. 

Consider the eigenvalue problem (jSJ)-©- Rescalling (jSJ) by the factor — , A = a 2 R 

A 

the eigenvalue problem can be written in the form Aw — XKw = 0, with 

A= { (D 2 -a 2 )J' K= \Nx-l 0j- (12) 
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Here w G T>a, with T>(A) the definition domain of the matricial differential operator 
A given by 

V A := [w = (W, 6) G (l 2 ( - i), i) V = £W = 6 = at « = -i -}. 

The following convergence result was proved. 

Theorem 1. JBj Let X be a parameter in the equation 

Aw - XKw = 0, (13) 

where A and K are linear operators, and the domain of A, Da, is a linear manifold 
that is dense in a Hilbert space H with the inner product (•, •). Let Da be contained 
in the domain of K , and assume that the following conditions are fulfilled: 

a) the operator A is a positive-definite, self adjoint operator; that is (Au,u) > 
and (Au, v) — (v, Au) = 0; 

b) the operator A~ l K can be extended to a completely continuous operator on the 
Hilbert space H n , where H n is the completion of Da under the norm (Au, u) 1 ^ 2 . 

Then the Galerkin method for calculating the eigenvalues of ( Tl3j) is a convergent 
process in H n . 

Using the definitions of the matricial differential operators, all the conditions of 
the theorem are satisfied so, the Galerkin method for computing the eigenvalues of 
dHD converges in the norm of H , with H the Hilbert space obtained by completing 
T>a (which is a preHilbert space) to a Hilbert space. 

Remark. Similarly, the convergence can be proved in the case of L 2 (0, 1). 

3 Galerkin type spectral methods 

The expansion functions used for the unknown fields encountered in various con- 
vection problems from hydrodynamic stability theory must have a basic property: 
they must be easy to evaluate. Trigonometric and polynomial functions have this 
property. A second requirement is the completeness of the sets of expansion func- 
tions. This assures that each function of the given space can be written as a linear 
combination of functions from the considered set (or, more likely, as a limit of such 
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a linear combination). The Chebyshev polynomials, the Legendre polynomials, the 
Hermite functions, the sine and cosine functions, satisfy this condition. 

In the Galerkin approach used here the basis (trial) functions satisfy the bound- 
ary conditions. In this case, following Rama Rao [11], the simplest choice seems to 
be to write W and as 

oo oo 

W = 2j a m h lm (x), 6 = 2J b m h 2 m(x), (14) 

m=0 m=0 

where h\ m {x) = (1 — Ax 2 ) m+2 , h 2m (x) = (1 — Ax 2 ) m+1 . With this choice, the un- 
known functions W and satisfy the boundary conditions (j7j). Replacing these 
expressions in and imposing the condition that the obtained equations are or- 
thogonal to hi n (x) and h 2n (x) respectively, n G N we obtained an algebraic system 
in the unknown coefficients a m and b m . The condition that these coefficients are 
nonnull gives us the secular (dispersion) equation. However, an important remark is 
in order: the physical parameter N representing the heating (cooling) rate is missing 
from this equation. 

Let us mention that the method is used in Ramma Rao [TT] in a convective 
instability problem of a heat conducting micropolar fluid layer situated between two 
rigid boundaries. In order to investigate the critical values of the Rayleigh number 
at which instability sets in the most rough approximation is taken, with only one 
term for each expression from (JHJ), so the approximate values of R are also crude. 
Nevertheless, in our case, for this approximation, in the classical case of Benard 
convection, corresponding to N = 0, the critical value of the Rayleigh number R is 
R = 1705.715 for a = 3.17, which is a very good approximation compared to the 
well-known value from Chandrasekhar [2j . We can conclude that this approximation 
works with good results only in the classical case. 

A mathematical explanation for the absence of the parameter N could be that the 
chosen set of expansion functions introduced an extraparity in the problem, leading 
to the loss of one of the physical parameter, in this case the heating (cooling) rate 
N. 

In [5] we considered also a basis of some hyperbolic functions for the expansion 

00 

of the unknown function W, i.e. W = W^CnW' For this choice the physical 

n=l 

parameter N was not present in the dispersion equation. This is why, we assume 
that a more suitable choice is to consider the general case 

00 

W{x) = Y,W 1 n C n {x) + W 2 n S n {x) 

n=l 
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where C n and S n are the Chandrasekhar sets of functions defined in [2] 

cosh(A n z) cos(\ n z) 



cosh(A n /2) cos(A n /2)' 

sinh(/i n z) sin(// n z) 
sinh(/i n /2) sin(/i n /2)' 



(15) 



with A n , /i n given in [2J by explicit values for n = 1,2,3,4 and by a recurrence 
relation for n > 4. 

From f6} 2 we obtain the expression of the unknown function 0, 



S(x) 



l R 0j(x) + A cosh(ax) + B sinh(ax) 



with A,B deduced from the boundary conditions q( ± =0. The functions 
Bi(x), % = 1,2, ...,4, depending on the coefficients and W 2 , have the form 

n = y. / ^(iVx- l)cosh(A w g) _ 2ATA n W„ 1 sinh(A ra x) 
n til (A 2 -a 2 )cosh(A n /2) (A 2 - a 2 ) 2 cosh(A n /2) 



© 2 (x) = E { 



W*(iVz - 1) cos(A„x) _ 2N\ n W^ sin(A n x) i 
\\ (A 2 + a 2 ) cos(A n /2) ~ (A 2 + a 2 ) 2 cos(A„/2) J ' 

^ / ^ 2 (iVx - 1) sinh(/x n g) _ 2jV/i ra iy TO 2 cosh(/x n x) -i 
" 3W ibl (/x 2 - a 2 ) sinh(/i n /2) 02)2^(^/2) J ' 



9 4 (x) = £ 



j W 2 (iVx — 1) sin(/i n x) 2iV/i n W / ' 2 cos(/i„ 



j; 



„=i I (/x 2 + a 2 ) sin(/x n /2) (/x 2 + a 2 ) 2 sin(// n /2) 

Let us replace this expression in The orthogonality relation on C m , S m , 

m G N imposed by the Galerkin procedure led us to an algebraic system for the 
unknown coefficients Wl and W 2 



oo 4 2 

E Wn[(An + « 4 )5n m - 2a 2 T nm ] - 2a 2 W 2 U nm = E Ce, + E Cm5 

n=l i=l fc=l 

oo 4 2 



(16) 



E -2a 2 V nm W l n + + a 4 )5 nm ] - 2a 2 P nm = E S 0l + E 3 



k 
n • 



n=l 



fc=l 



with 



T 



and 



r 1 



5-1 



1/2 



1/2 
1/2 



c, 



e, 



1/2 
1/2 



cosh(ax)C m (x); C 5 
cosh(ax)5' m (x); S 2 



1/2 



-1/2 
1/2 



6i(x)C m (a;); 5*6, 



-1/2 



1/2 



1/2 



sinh(ax)C m (x) 
sinh(aa;)S' m (a;); 



Qi(x)S m (x) 



-1/2 



This time, the secular equation depends on N and it follows from the condition that 
not all these coefficients vanish. Numerical values of the Rayleigh number are then 
obtained and displayed in Table 1 in comparison with previous results. 



N 


a 2 


Ra-m 


R a — here 





9.711 


1715.079324 


1708.54 


1 


9.711 


1711.742588 


1651.04 


2 


9.711 


1701.891001 


1609.12 


1 


10.0 


1712.257687 


1651.1 


4 


10.0 


1664.341789 


1560.8 


4 


12.0 


1685.422373 


1739.2 


10 


9.0 


1482.527042 


1366.02 


11 


9.0 


1446.915467 


1366.05 


12 


9.00 


1411.401914 


1354.7 



Table 1. Numerical evaluations of the Rayleigh number for various values of the 

parameters N and a. 

For the eigenvalue problem ([H])-®, in [5], in order to avoid the loss of the pa- 
rameter N different sets of orthogonal functions based on polynomials, namely on 
shifted Legendre polynomials (SLP) on [0, 1] were proposed. The method is similar 

to the one presented here. Instead of {hi m (x)} m and {h 2m (x)} m from L 2 (^ — -, 

we used the orthogonal sets from L 2 (0, 1), 



JO 



Q m+1 (t)dtds 



Qm+3 Q 



m+1 



Q 



m+l 



Q 



m—l 



(2m + 3) (2m + 5) (2m + 1) (2m + 3) . 



{<p m (z)} m : (j) m {z) = [ Q m (t)dt = %±J 

Jo 2 (2m + 1) 

respectively, with Q m the classical Legendre polynomials defined on [—1, 1]. 

In this case, the expression of the secular equation contains the physical param- 
eter N, so good numerical evaluations of the Rayleigh number for various values of 
N and a were obtained. 

In [6], a general Galerkin type method is proposed for the problem written in 
the general form (fTDI) -(lT]). The unknown function is written as a Fourier series [6] 
of the form 

oo 

@ = ^2 Am C0S iPmz) + B m sin(q m z), (17) 

m=l 

where p m = (2m — l)ir, q m = 1mi\ which implies that O satisfies the boundary 
conditions (02). The expression of O, introduced in (flOl) 1 leads to an expression of 

oo 

W in the form W = A m f m (z) + B m g m (z) in which the boundary conditions ([7]) 

m=l 

are also considered in order to find A m and B m . However, in our case, the function 
f(z) is a constant one and the application of the method in this form to ([6])-([7j) does 
not lead to a correct expression of W . 



4 Conclusion 

In this paper a problem of convection with uniform internal heat source is in- 
vestigated. We complete a previous analytical study [1], [5] with some comments 
on the choice of the expansion functions and their importance for the convergence 
of the Galerkin method. The importance of the form of the system of ordinary 
differential equations which describe the eigenvalue problem governing to the linear 
stability of the stationary solution with respect to this choice is pointed out. We 
present numerical results for the new introduced methods which are similar to the 
ones obtained before. 

The main conclusion of our analytical and numerical study performed in this 
paper and the previous ones is that the choice of subspaces of trial functions with 
respect to whom the approximation problems are solved influences the form of the 
algebraic system and also the numerical evaluations. The good numerical results 
obtained for small values of the spectral parameter are justified by the accuracy of 
spectral methods. 
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